function theeijk(sigma)
format long

hermithd{1}=@(x) 1+x-x;
hermithd{2}=@(x) x/sigma;
hermithd{3}=@(x) (x.^2/sigma^2-1)/sqrt(2);
hermithd{4}=@(x) (x.^3/sigma^3-3*x/sigma)/sqrt(6);
hermithd{5}=@(x) (x.^4/sigma^4-6*x.^2/sigma^2+3)/sqrt(24);

gaussian =@(x) 1/sqrt(2*pi)/sigma*exp(-x.^2/2/sigma^2);

testij=zeros(5);
hermitijk=zeros(5,5,5);

for ii=1:5
    for jj=1:5
        tmphd=@(x) hermithd{ii}(x).*hermithd{jj}(x).*gaussian(x);
        testij(ii,jj)=quadgk(tmphd,-10,10,'RelTol',1e-12,'AbsTol',1e-12);
    end
end

for ii=1:5
    for jj=1:5
        for kk=1:5
            tmphd =@(x) hermithd{ii}(x).*hermithd{jj}(x).*hermithd{kk}(x).*gaussian(x);
            hermitijk(ii,jj,kk)=quadgk(tmphd,-10,10,'RelTol',1e-12,'AbsTol',1e-12);
        end
    end
end

fname='hermitdata';
save(fname, 'hermitijk', 'testij', 'hermithd','gaussian');